MODULE SRFSN_LWIMPSAD_MOD
CONTAINS
SUBROUTINE SRFSN_LWIMPSAD(KIDIA ,KFDIA    ,KLON   ,LDLAND, &
 & PTMST    ,LDREGSF  ,&
 & PSSNM1M5 ,PTSNM1M5 ,PRSNM1M  ,PTSAM1M5 ,PTIAM1M5 ,PHLICEM1M5,&
 & PSLRFL5  ,PSSRFLTI5,PFRTI    ,PAHFSTI5 ,PEVAPTI5 ,&
 & PSSFC5   ,PSSFL5   ,PEVAPSNW5,&
 & PTSFC5   ,PTSFL5   ,&
 & YDCST    ,YDVEG    ,YDSOIL   ,YDFLAKE  , YDURB,&
 & PTSN5    ,PGSN5    ,&                             
 & PSSNM1M  ,PTSNM1M  ,PTSAM1M  ,PTIAM1M,&
 & PSLRFL   ,PSSRFLTI ,PAHFSTI  ,PEVAPTI,&
 & PSSFC    ,PSSFL    ,PEVAPSNW ,&
 & PTSFC    ,PTSFL    ,&
 & PTSN    , PGSN )

USE PARKIND1  , ONLY : JPIM, JPRB
USE YOMHOOK   , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_CST   , ONLY : TCST
USE YOS_VEG   , ONLY : TVEG
USE YOS_SOIL  , ONLY : TSOIL
USE YOS_FLAKE , ONLY : TFLAKE
USE YOS_URB   , ONLY : TURB

#ifdef DOC
! (C) Copyright 2011- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!**** *SRFSNS_LWIMPAD* - CONTAINS SNOW PARAMETRIZATION
!                        (Adjoint)
!
!     PURPOSE.
!     --------
!          COMPUTES CHANGES IN SNOW TEMPERATURE

!**   INTERFACE.
!     ----------
!          *SRFSN_SLWIMPAD* IS CALLED FROM *SURFTSTPSAD*.

!     PARAMETER   DESCRIPTION                                          UNITS
!     ---------   -----------                                          -----

!     INPUT PARAMETERS (INTEGER):
!     *KIDIA*      START POINT
!     *KFDIA*      END POINT
!     *KLON*       NUMBER OF GRID POINTS PER PACKET

!     INPUT PARAMETERS (REAL):
!     *PTMST*      TIME STEP                                            S
!     *PFRTI*      TILE FRACTIONS                                       -
 
!     INPUT PARAMETERS (LOGICAL):
!     *LDREGSF*    TRUE WHEN REGULARIZATION USED

!     INPUT PARAMETERS AT T-1 OR CONSTANT IN TIME (REAL):
!  Trajectory  Perturbation  Description                                Unit
!  PSSNM1M5    PSSNM1M       SNOW MASS (per unit area)                  kg/m2/s
!  PTSNM1M5    PTSNM1M       SNOW TEMPERATURE                           K
!  PRSNM1M     ---           SNOW DENSITY                               kg/m3 
!  PTSAM1M5    PTSAM1M       SOIL TEMPERATURE                           K
!  PHLICEM1M5  ---           LAKE ICE THICKNESS                         m
!  PSLRFL5     PSLRFL        NET LONGWAVE  RADIATION AT THE SURFACE     W/m2
!  PSSRFLTI5   PSSRFLTI      NET SHORTWAVE RADIATION AT THE SURFACE,    W/m2
!                            FOR EACH TILE  
!  PAHFSTI5    PAHFSTI       TILE SENSIBLE HEAT FLUX                    W/m2
!  PEVAPTI5    PEVAPTI       TILE EVAPORATION                           kg/m2/s
!  PSSFC5      PSSFC         CONVECTIVE  SNOW FLUX AT THE SURFACE       kg/m2/s
!  PSSFL5      PSSFL         LARGE SCALE SNOW FLUX AT THE SURFACE       kg/m2/s
!  PEVAPSNW5   PEVAPSNW      EVAPORATION FROM SNOW UNDER FOREST         kg/m2/s
!  PTSFC5      PTSFC         Convective Throughfall at the surface      kg/m2/s
!  PTSFL5      PTSFL         Large Scale Throughfall at the surface     kg/m2/s

!     PARAMETERS AT T+1 :
!  Trajectory  Perturbation  Description                                Unit
!  PTSN5       PTSN          SNOW TEMPERATURE                           K

!    FLUXES FROM SNOW SCHEME:
!  Trajectory  Perturbation  Description                                Unit
!  PGSN5       PGSN          GROUND HEAT FLUX FROM SNOW DECK TO SOIL    W/m2 (#)


! (#) THOSE TWO QUANTITIES REPRESENT THE WHOLE GRID-BOX. IN RELATION
!       TO THE DOCUMENTATION, THEY ARE PGSN=Fr_s*G_s

!     METHOD.
!     -------
!     Based on the original snow (as in ERA-40) with the following updates:
!     - Liquid water as a diagnostics
!     - Interception of rainfall

!     EXTERNALS.
!     ----------

!     REFERENCE.
!     ----------
!          SEE SOIL PROCESSES' PART OF THE MODEL'S DOCUMENTATION FOR
!     DETAILS ABOUT THE MATHEMATICS OF THIS ROUTINE.

!     Original:
!     ---------
!     M. Janiskova              E.C.M.W.F.     28-03-2012  

!     Modifications
!     -------------

!     ------------------------------------------------------------------
#endif

IMPLICIT NONE

! Declaration of arguments

INTEGER(KIND=JPIM), INTENT(IN)   :: KIDIA
INTEGER(KIND=JPIM), INTENT(IN)   :: KFDIA
INTEGER(KIND=JPIM), INTENT(IN)   :: KLON

LOGICAL,            INTENT(IN)   :: LDREGSF
LOGICAL           ,INTENT(IN)    :: LDLAND(:)

REAL(KIND=JPRB),    INTENT(IN)   :: PTMST
REAL(KIND=JPRB),    INTENT(IN)   :: PSSNM1M5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PTSNM1M5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PRSNM1M(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PTSAM1M5(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PTIAM1M5(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PHLICEM1M5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSLRFL5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSRFLTI5(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PFRTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PAHFSTI5(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PEVAPTI5(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSFC5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSFL5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PEVAPSNW5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PTSFC5(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PTSFL5(:)
TYPE(TCST),         INTENT(IN)   :: YDCST
TYPE(TVEG),         INTENT(IN)   :: YDVEG
TYPE(TSOIL),        INTENT(IN)   :: YDSOIL
TYPE(TFLAKE),       INTENT(IN)   :: YDFLAKE
TYPE(TURB),         INTENT(IN)   :: YDURB

REAL(KIND=JPRB),    INTENT(OUT)  :: PTSN5(:)
REAL(KIND=JPRB),    INTENT(OUT)  :: PGSN5(:)

REAL(KIND=JPRB),    INTENT(INOUT) :: PSSNM1M(:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PTSNM1M(:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PTSAM1M(:,:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PTIAM1M(:,:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PSLRFL(:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PSSRFLTI(:,:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PAHFSTI(:,:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PEVAPTI(:,:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PSSFC(:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PSSFL(:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PEVAPSNW(:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PTSFC(:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PTSFL(:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PTSN(:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PGSN(:)


!      LOCAL VARIABLES

LOGICAL            :: LLNOSNOW(KLON)

REAL(KIND=JPRB),DIMENSION(KLON) :: ZSSFC5,ZSSFL5
REAL(KIND=JPRB),DIMENSION(KLON) :: ZTSFC5,ZTSFL5,ZLWC5,ZDSNR5,ZDSNRP5
REAL(KIND=JPRB),DIMENSION(KLON) :: ZSSTAR5,ZSSTARC5,ZHFLUX5,ZLWC15
REAL(KIND=JPRB),DIMENSION(KLON) :: ZFUNC5,ZDFUNC5,ZCONSA5,ZCONSB5,ZCONS25
REAL(KIND=JPRB),DIMENSION(KLON) :: ZCOND15,ZCOND25,ZCOND35,ZCOND45
REAL(KIND=JPRB),DIMENSION(KLON) :: ZCOND55,ZCOND65
REAL(KIND=JPRB),DIMENSION(KLON) :: ZRS5,ZHFLUXRF15,ZHFLUXRF25,ZSSTAR15
REAL(KIND=JPRB),DIMENSION(KLON) :: ZTSTARI5,ZTSTAR5,ZDTH5,ZDTHI5,ZDTM5,ZMSN5
REAL(KIND=JPRB) :: ZHFLUXP5
REAL(KIND=JPRB) :: ZSNRES5,ZWFLUXRF5
REAL(KIND=JPRB) :: ZGSN5,ZDSN5

REAL(KIND=JPRB),DIMENSION(KLON) :: ZLICE,ZFRSN,ZIFRSN
REAL(KIND=JPRB),DIMENSION(KLON) :: ZSSFC,ZSSFL,ZLWC,ZDSNR,ZDSNRP,ZTSFC,ZTSFL
REAL(KIND=JPRB),DIMENSION(KLON) :: ZGRIDFRAC,ZCONSLWC
REAL(KIND=JPRB),DIMENSION(KLON) :: ZCONSAMAX
REAL(KIND=JPRB) :: ZFUNC,ZDFUNC,ZDSN,ZDTH,ZDTHI,ZDTM,ZEXPF,&
 &        ZGSN,ZHFLUXP,ZHFLUX,ZLINA,ZMSN,ZRS,ZSNRES,ZSOILRES,ZSSTAR,ZSSTAR1,&
 &        ZT0,ZTSTAR,ZTSTARI,ZHOICE,ZTMST,ZCONSA,ZCONSB,ZCONS2,&
 &        ZWFLUXRF,ZHFLUXRF1,ZITEMPAMP
REAL(KIND=JPRB) :: ZCOND1,ZCOND2,ZCOND3,ZCOND4,ZCOND5,ZCOND6
REAL(KIND=JPRB),DIMENSION(KLON) :: ZTSURF
REAL(KIND=JPRB),DIMENSION(KLON,1) :: ZTSURF5

REAL(KIND=JPRB) :: ZRDSNMAX, ZRDSNRESMAX

INTEGER(KIND=JPIM) :: JL

REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!     ------------------------------------------------------------------

IF (LHOOK) CALL DR_HOOK('SRFSN_LWIMPSAD_MOD:SRFSN_LWIMPSAD',0,ZHOOK_HANDLE)
ASSOCIATE(RDAY=>YDCST%RDAY, RLMLT=>YDCST%RLMLT, RLSTT=>YDCST%RLSTT, &
 & RLVTT=>YDCST%RLVTT, RPI=>YDCST%RPI, RTT=>YDCST%RTT, &
 & LEFLAKE=>YDFLAKE%LEFLAKE, LEURBAN=>YDURB%LEURBAN, RH_ICE_MIN_FLK=>YDFLAKE%RH_ICE_MIN_FLK, &
 & RALAMSN=>YDSOIL%RALAMSN, RDSNMAX=>YDSOIL%RDSNMAX, RFRSMALL=>YDSOIL%RFRSMALL, &
 & RFRTINY=>YDSOIL%RFRTINY, RHOCI=>YDSOIL%RHOCI, RHOICE=>YDSOIL%RHOICE, &
 & RLAMICE=>YDSOIL%RLAMICE, RLWCSWEA=>YDSOIL%RLWCSWEA, &
 & RLWCSWEB=>YDSOIL%RLWCSWEB, RLWCSWEC=>YDSOIL%RLWCSWEC, RTAUA=>YDSOIL%RTAUA, &
 & RTAUF=>YDSOIL%RTAUF, RTEMPAMP=>YDSOIL%RTEMPAMP, &
 & RVLAMSK_DESERT=>YDVEG%RVLAMSK_DESERT, RVLAMSK_SNOW=>YDVEG%RVLAMSK_SNOW)

!*         1.    SET UP SOME CONSTANTS.
!                --- -- ---- ----------
!*               PHYSICAL CONSTANTS.
!                -------- ----------
! RESISTANCE FOR HEAT CONDUCTION IN HALF TOP SOIL LAYER
!  (ESTIMATE FROM SKIN LAYER CONDUCTIVITY FOR BARE SOIL; W/M2K)


ZHOICE=1.0_JPRB/RHOICE
ZT0=RTT
ZTMST=1.0_JPRB/PTMST
ZEXPF=EXP(-RTAUF*PTMST/RDAY)
ZLINA=RTAUA*PTMST/RDAY
ZITEMPAMP=1.0_JPRB/RTEMPAMP


ZRDSNMAX=1.0_JPRB
ZRDSNRESMAX=1._JPRB

!*   FIND LAKE POINTS WITH ICE COVER
DO JL=KIDIA,KFDIA
  IF ( PHLICEM1M5(JL) > RH_ICE_MIN_FLK  .AND. LEFLAKE ) THEN
    ZLICE(JL)=1._JPRB
  ELSE
    ZLICE(JL)=0._JPRB
  ENDIF
ENDDO

!     ------------------------------------------------------------------
!*         2. NEW SNOW T AND MASS INCLUDING GROUND HEAT FLUX AND MELTING.
!             -----------------------------------------------------------

DO JL=KIDIA,KFDIA
  ZFRSN(JL)=MAX(PFRTI(JL,5)+PFRTI(JL,7),RFRTINY)
  ZIFRSN(JL)=1.0_JPRB/ZFRSN(JL)
  IF (ZFRSN(JL) < RFRSMALL) THEN
    LLNOSNOW(JL)=.TRUE.
  ELSE
    LLNOSNOW(JL)=.FALSE.
  ENDIF
  ZGRIDFRAC(JL) = (PFRTI(JL,3)+PFRTI(JL,4)+PFRTI(JL,5)&
   & +PFRTI(JL,6)+PFRTI(JL,7)+PFRTI(JL,8))

  IF ( LEURBAN ) THEN
   ZGRIDFRAC=(PFRTI(JL,3)+PFRTI(JL,4)+PFRTI(JL,5)&
   & +PFRTI(JL,6)+PFRTI(JL,7)+PFRTI(JL,8)+PFRTI(JL,10))
  ENDIF


  IF ( LEFLAKE ) THEN
    IF ( PFRTI(JL,9) .EQ. 1._JPRB ) THEN
      ZGRIDFRAC(JL) = 0._JPRB ! FULL COVER OF LAKE --NO SNOW THERE !
    ELSE
      ZGRIDFRAC(JL) = ZGRIDFRAC(JL)+PFRTI(JL,9)
    ENDIF
  ENDIF

  ZSSFC5(JL) = ZGRIDFRAC(JL)*PSSFC5(JL)
  ZSSFL5(JL) = ZGRIDFRAC(JL)*PSSFL5(JL)
  ZTSFC5(JL) = ZGRIDFRAC(JL)*PTSFC5(JL)
  ZTSFL5(JL) = ZGRIDFRAC(JL)*PTSFL5(JL)
ENDDO

DO JL=KIDIA,KFDIA
  IF (LLNOSNOW(JL)) THEN
    ZSSTAR5(JL) = PSSNM1M5(JL)+PTMST*(ZSSFC5(JL)+ZSSFL5(JL))
    PTSN5(JL) = ZT0
    PGSN5(JL) = 0.0_JPRB
    ZLWC5(JL) = 0.0_JPRB
  ELSE
    ! Snow on land or ice: select bottom boundary condition
    IF (LDLAND(JL)) THEN
      ZSOILRES=1.0_JPRB/RVLAMSK_DESERT
      ZTSURF5(JL,1)=PTSAM1M5(JL,1)
    ELSE
      IF(YDSOIL%LESNICE)THEN
        ZSOILRES=1.0_JPRB/RVLAMSK_SNOW
        ZTSURF5(JL,1)=PTIAM1M5(JL,1)
      ELSE
        ZSOILRES=1.0_JPRB/RVLAMSK_DESERT
        ZTSURF5(JL,1)=PTSAM1M5(JL,1)
      ENDIF
    ENDIF

!           NET HEAT FLUX AT SNOW SURFACE; THIS IS MULTIPLIED
!           BY FRACTION BECAUSE EQUATIONS APPLY TO TOTAL SNOW MASS
!           IN THE GRID SQUARE

    ZHFLUXP5 = PFRTI(JL,5)*(PAHFSTI5(JL,5)+RLSTT*PEVAPTI5(JL,5))&
     & + PFRTI(JL,7)*(PAHFSTI5(JL,7)+RLVTT*(PEVAPTI5(JL,7)-PEVAPSNW5(JL))&
     & + RLSTT*PEVAPSNW5(JL))&
     & + PFRTI(JL,5)*PSSRFLTI5(JL,5)&
     & + PFRTI(JL,7)*PSSRFLTI5(JL,7)&
     & + (PFRTI(JL,5)+PFRTI(JL,7))*PSLRFL5(JL)
    ZHFLUX5(JL) = ZIFRSN(JL)*ZHFLUXP5

! PRELIMINARY SNOW MASS
    ZSSTAR5(JL) = PSSNM1M5(JL)+PTMST*&
     & (ZSSFC5(JL)+ZSSFL5(JL)+PFRTI(JL,5)*PEVAPTI5(JL,5)&
     & +PFRTI(JL,7)*PEVAPSNW5(JL))
    ! CORRECT FOR NEGATIVE VALUES OF SNOW MASS AFTER P-E
    ZSSTARC5(JL) = ZSSTAR5(JL)
    IF (ZSSTARC5(JL) < 0.0_JPRB) THEN
!      ZSSTAR=MAX(ZSSTAR,0.0_JPRB)
      ZSSTAR5(JL) = 0.0_JPRB
    ENDIF

! SNOW LIQUID WATER - IMPLICIT
    ! SNOW DEPTH (REAL)
    ZDSNR5(JL) = PSSNM1M5(JL)/(PRSNM1M(JL)*ZFRSN(JL))

    ! SNOW LIQUID WATER CAPACITY
    ZCONSLWC(JL) = MAX(0._JPRB,(RLWCSWEA-PRSNM1M(JL))/RLWCSWEA)
    ZLWC5(JL) = PSSNM1M5(JL)*(RLWCSWEB+(RLWCSWEC-RLWCSWEB)*ZCONSLWC(JL))
    ZLWC15(JL) = ZLWC5(JL)

    ! LWC IS RESCALLED CONSIDERING EFFECTIVE DEPTH (RDSNMAX)
    ! ZLWC(JL)=ZLWC(JL)*MIN(ZDSNR(JL),ZRDSNMAX)/ZDSNR(JL)
    IF (ZDSNR5(JL) < ZRDSNMAX) THEN
      ZDSNRP5(JL) = ZDSNR5(JL)
    ELSE
      ZDSNRP5(JL) = ZRDSNMAX
    ENDIF
    ZLWC5(JL) = ZLWC15(JL)*ZDSNRP5(JL)/ZDSNR5(JL)  

    ! ANALYTICAL FUNCTIONS - SNOW LIQUID WATER
    IF ( PTSNM1M5(JL) < (ZT0-0.5_JPRB*RTEMPAMP) ) THEN
      ZFUNC5(JL) = 0._JPRB
      ZDFUNC5(JL) = 0._JPRB
    ELSEIF(PTSNM1M5(JL) < ZT0) THEN
      ZFUNC5(JL) = 1._JPRB+SIN(RPI*(PTSNM1M5(JL)-ZT0)*ZITEMPAMP)
      ZDFUNC5(JL) = COS(RPI*(PTSNM1M5(JL)-ZT0)*ZITEMPAMP)*RPI*ZITEMPAMP
    ELSE
      ZFUNC5(JL) = 0._JPRB
      ZDFUNC5(JL) = RPI*ZITEMPAMP
    ENDIF

! TERMS FOR IMPLICIT TEMPERATURE EQUATION
    ZCONSB5(JL) = RLMLT*ZLWC5(JL)*ZDFUNC5(JL)
    ZCONSAMAX(JL) =(PRSNM1M(JL)*RHOCI*ZRDSNMAX)*ZHOICE
    ZCOND15(JL) = (PSSNM1M5(JL)*RHOCI*ZHOICE)*ZIFRSN(JL)
    IF (ZCOND15(JL) < ZCONSAMAX(JL)) THEN
      ZCONSA5(JL) = ZCOND15(JL)
    ELSE
      ZCONSA5(JL) = ZCONSAMAX(JL)
    ENDIF
    ZCONS25(JL) = PTMST/(ZCONSA5(JL)+ZCONSB5(JL))

! REAL SNOW DEPTH
    ZCOND25(JL) = PSSNM1M5(JL)/(PRSNM1M(JL)*ZFRSN(JL))
    IF (ZCOND25(JL) < ZRDSNRESMAX) THEN
      ZDSN5 = ZCOND25(JL)
    ELSE
      ZDSN5 = ZRDSNRESMAX
    ! ZDSN5 = ZCOND25(JL)
    ENDIF
    ! RESISTANCE FOR HEAT FLUX BETWEEN SNOW AND SOIL LAYER
    ZSNRES5  = 0.5_JPRB*ZDSN5/(RLAMICE*(PRSNM1M(JL)*ZHOICE)**RALAMSN)
    ZRS5(JL) = 1.0_JPRB/(ZSNRES5+ZSOILRES)

!  ACCOUNT FOR THROUGHFALL INTERCEPTION
    ! mass advected
    ZWFLUXRF5 = (PFRTI(JL,5)+PFRTI(JL,7))*(ZTSFC5(JL)+ZTSFL5(JL))
    ! phase changes
    ZCOND35(JL) = ZWFLUXRF5*PTMST
    ZCOND45(JL) = ZLWC5(JL)*(1.0_JPRB-ZFUNC5(JL))
    IF (ZCOND35(JL) < ZCOND45(JL)) THEN
      ZWFLUXRF5 = ZCOND35(JL)
    ELSE
      ZWFLUXRF5 = ZCOND45(JL)
    ENDIF
    ZCOND55(JL) = ZCONSA5(JL)*(ZT0-PTSNM1M5(JL))*ZTMST+ZT0*ZRS5(JL)
    ZCOND65(JL) = RLMLT*ZWFLUXRF5*ZTMST
    IF (ZCOND55(JL) < ZCOND65(JL)) THEN
      ZHFLUXRF15(JL) = ZCOND55(JL)
    ELSE
      ZHFLUXRF15(JL) = ZCOND65(JL)
    ENDIF
    ZHFLUXRF25(JL) = ZHFLUXRF15(JL)
    IF (ZHFLUXRF25(JL) < 0.0_JPRB) THEN
      ZHFLUXRF15(JL) = 0.0_JPRB
    ENDIF
    !mass intercepted - recalculate from ZHFLUXRF1
    ZWFLUXRF5 = ZHFLUXRF15(JL)/RLMLT
    ZSSTAR5(JL) = ZSSTAR5(JL)+ZWFLUXRF5*PTMST
! SOLVE IMPLICIT TEMPERATURE EQUATION
    ZTSTARI5(JL) = 1.0_JPRB/(1.0_JPRB+ZCONS25(JL)*ZRS5(JL))
    ZTSTAR5(JL) = (PTSNM1M5(JL)+ZCONS25(JL) &
     & *(ZHFLUX5(JL)+ZHFLUXRF15(JL)+ZTSURF5(JL,1)*ZRS5(JL)))*ZTSTARI5(JL)
    IF (ZTSTAR5(JL) <= ZT0) THEN
! NO MELTING OF SNOW
      PTSN5(JL) = ZTSTAR5(JL)
      PGSN5(JL) = (ZTSTAR5(JL)-ZTSURF5(JL,1))*ZRS5(JL)
    ELSE
! MELTING OF SNOW:
! HEATING PART OF THE TIME STEP
      ZDTHI5(JL) = 1.0_JPRB/(ZHFLUX5(JL)+ZHFLUXRF15(JL) &
       & -(ZT0-ZTSURF5(JL,1))*ZRS5(JL))
      ZDTH5(JL) = ((ZT0-PTSNM1M5(JL))*(ZCONSA5(JL)+ZCONSB5(JL)))*ZDTHI5(JL)
! MELTING PART OF THE TIME STEP
      ZDTM5(JL) = PTMST-ZDTH5(JL)
      ZGSN5 = (ZT0-ZTSURF5(JL,1))*ZRS5(JL)
      ZMSN5(JL) = (ZHFLUX5(JL)+ZHFLUXRF15(JL)-ZGSN5)/RLMLT
      ZSSTAR15(JL) = ZSSTAR5(JL)-ZDTM5(JL)*ZMSN5(JL)*ZFRSN(JL)
      IF (ZSSTAR15(JL) >= 0.0_JPRB) THEN
! NOT ALL THE SNOW HAS BEEN MELTED
        PTSN5(JL) = ZT0
        PGSN5(JL) = (ZT0-ZTSURF5(JL,1))*ZRS5(JL)
      ELSE
! ALL THE SNOW HAS BEEN MELTED; COMPUTE TIME IT TOOK TO MELT
        PTSN5(JL) = ZT0
        PGSN5(JL) = ( ((ZT0-ZTSURF5(JL,1))*ZRS5(JL))*(ZDTM5(JL)+ZDTH5(JL))&
         & +(ZHFLUX5(JL)+ZHFLUXRF15(JL))*(PTMST-ZDTM5(JL)-ZDTH5(JL)) )*ZTMST
      ENDIF
    ENDIF

  ENDIF
ENDDO

!*         5. NORMALIZE QUANTITIES TO THE GRID-SQUARE.
!             ----------------------------------------

DO JL=KIDIA,KFDIA
  PGSN5(JL) = ZFRSN(JL)*PGSN5(JL)
ENDDO

!           0.  ADJOINT CALCULATIONS
!               --------------------

!* Set local variables to zero

ZSSFC(:)  = 0.0_JPRB
ZSSFL(:)  = 0.0_JPRB
ZLWC(:)   = 0.0_JPRB
ZDSNR(:)  = 0.0_JPRB
ZDSNRP(:) = 0.0_JPRB
ZTSFC(:)  = 0.0_JPRB
ZTSFL(:)  = 0.0_JPRB

!*         0. 5. NORMALIZE QUANTITIES TO THE GRID-SQUARE.
!                ----------------------------------------

DO JL=KIDIA,KFDIA
  PGSN(JL) = ZFRSN(JL)*PGSN (JL)
ENDDO

!     ------------------------------------------------------------------
!*         0. 2. NEW SNOW T AND MASS INCLUDING GROUND HEAT FLUX AND MELTING.
!                -----------------------------------------------------------

DO JL=KIDIA,KFDIA
  IF (LLNOSNOW(JL)) THEN
    ZSSTAR = 0.0_JPRB

    ZLWC(JL) = 0.0_JPRB
    PGSN(JL) = 0.0_JPRB
    PTSN(JL) = 0.0_JPRB
    PSSNM1M(JL) = PSSNM1M(JL)+ZSSTAR
    ZSSFC(JL) = ZSSFC(JL)+PTMST*ZSSTAR
    ZSSFL(JL) = ZSSFL(JL)+PTMST*ZSSTAR
  ELSE
    ZTSTAR = 0.0_JPRB
    ZRS  = 0.0_JPRB
    ZDTM = 0.0_JPRB
    ZDTH = 0.0_JPRB
    ZHFLUX = 0.0_JPRB
    ZHFLUXRF1 = 0.0_JPRB
    ZMSN = 0.0_JPRB
    ZGSN = 0.0_JPRB
    ZSSTAR1 = 0.0_JPRB
    ZDTHI  = 0.0_JPRB
    ZSSTAR = 0.0_JPRB
    ZDTHI  = 0.0_JPRB
    ZTSTARI = 0.0_JPRB
    ZCONS2 = 0.0_JPRB
    ZWFLUXRF = 0.0_JPRB
    ZCOND5 = 0.0_JPRB
    ZCOND6 = 0.0_JPRB
    ZCONSA = 0.0_JPRB
    ZCOND3 = 0.0_JPRB
    ZCOND4 = 0.0_JPRB
    ZFUNC  = 0.0_JPRB
    ZSNRES = 0.0_JPRB
    ZDSN   = 0.0_JPRB
    ZCOND2 = 0.0_JPRB
    ZCOND1 = 0.0_JPRB
    ZCONS2 = 0.0_JPRB
    ZCONSA = 0.0_JPRB
    ZCONSB = 0.0_JPRB
    ZDFUNC = 0.0_JPRB
    ZHFLUXP = 0.0_JPRB
    ZTSURF(JL)=0.0_JPRB
    IF (ZTSTAR5(JL) <= ZT0) THEN
! NO MELTING OF SNOW
      ZRS = ZRS+(ZTSTAR5(JL)-ZTSURF5(JL,1))*PGSN(JL)
      ZTSTAR = ZTSTAR+ZRS5(JL)*PGSN(JL)
      !*ZTSURF(JL) = ZTSURF(JL)-ZRS5(JL)*PGSN(JL)
      IF(YDSOIL%LESNICE)THEN
        PTIAM1M(JL,1) = PTIAM1M(JL,1)-ZRS5(JL)*PGSN(JL)
      ELSE
        PTSAM1M(JL,1) = PTSAM1M(JL,1)-ZRS5(JL)*PGSN(JL)
      ENDIF
      PGSN(JL) = 0.0_JPRB
      ZTSTAR = ZTSTAR+PTSN(JL)
      PTSN(JL) = 0.0_JPRB
    ELSE
! MELTING OF SNOW:
      IF (ZSSTAR15(JL) >= 0.0_JPRB) THEN
! NOT ALL THE SNOW HAS BEEN MELTED
        ZRS = ZRS+(ZT0-ZTSURF5(JL,1))*PGSN(JL)
        !*ZTSURF(JL) = ZTSURF(JL)-ZRS5(JL)*PGSN(JL)
        IF(YDSOIL%LESNICE)THEN
          PTIAM1M(JL,1) = PTIAM1M(JL,1)-ZRS5(JL)*PGSN(JL)
        ELSE
          PTSAM1M(JL,1) = PTSAM1M(JL,1)-ZRS5(JL)*PGSN(JL)
        ENDIF
        PGSN(JL) = 0.0_JPRB
        PTSN(JL) = 0.0_JPRB
      ELSE
! ALL THE SNOW HAS BEEN MELTED; COMPUTE TIME IT TOOK TO MELT
        ZDTM = ZDTM+ZTMST*((ZT0-ZTSURF5(JL,1))*ZRS5(JL))*PGSN(JL)
        ZDTH = ZDTH+ZTMST*((ZT0-ZTSURF5(JL,1))*ZRS5(JL))*PGSN(JL)
        ZRS = ZRS+ZTMST*(ZT0-ZTSURF5(JL,1))*(ZDTM5(JL)+ZDTH5(JL))*PGSN(JL)
        !*ZTSURF(JL) = ZTSURF(JL) &
        !* & -ZTMST*ZRS5(JL)*(ZDTM5(JL)+ZDTH5(JL))*PGSN(JL)
        IF(YDSOIL%LESNICE)THEN
          PTIAM1M(JL,1) = PTIAM1M(JL,1) &
           & -ZTMST*ZRS5(JL)*(ZDTM5(JL)+ZDTH5(JL))*PGSN(JL)
        ELSE
          PTSAM1M(JL,1) = PTSAM1M(JL,1) &
           & -ZTMST*ZRS5(JL)*(ZDTM5(JL)+ZDTH5(JL))*PGSN(JL)
        ENDIF
        ZDTM = ZDTM-ZTMST*(ZHFLUX5(JL)+ZHFLUXRF15(JL))*PGSN(JL)
        ZDTH = ZDTH-ZTMST*(ZHFLUX5(JL)+ZHFLUXRF15(JL))*PGSN(JL)
        ZHFLUX = ZHFLUX+ZTMST*(PTMST-ZDTM5(JL)-ZDTH5(JL))*PGSN(JL)
        ZHFLUXRF1 = ZHFLUXRF1+ZTMST*(PTMST-ZDTM5(JL)-ZDTH5(JL))*PGSN(JL)
        PGSN(JL) = 0.0_JPRB
        PTSN(JL) = 0.0_JPRB
      ENDIF
! MELTING PART OF THE TIME STEP
      ZSSTAR = ZSSTAR+ZSSTAR1
      ZMSN = ZMSN-ZFRSN(JL)*ZDTM5(JL)*ZSSTAR1
      ZDTM = ZDTM+ZFRSN(JL)*ZMSN5(JL)*ZSSTAR1

      ZHFLUX = ZHFLUX+ZMSN/RLMLT
      ZHFLUXRF1 = ZHFLUXRF1+ZMSN/RLMLT
      ZGSN = ZGSN-ZMSN/RLMLT

      ZRS = ZRS+(ZT0-ZTSURF5(JL,1))*ZGSN
      !*ZTSURF(JL) = ZTSURF(JL)-ZRS5(JL)*ZGSN
      IF(YDSOIL%LESNICE)THEN
        PTIAM1M(JL,1) = PTIAM1M(JL,1)-ZRS5(JL)*ZGSN
      ELSE
        PTSAM1M(JL,1) = PTSAM1M(JL,1)-ZRS5(JL)*ZGSN
      ENDIF
      ZDTH = ZDTH-ZDTM
! HEATING PART OF THE TIME STEP
      ZDTHI = ZDTHI+((ZT0-PTSNM1M5(JL))*(ZCONSA5(JL)+ZCONSB5(JL)))*ZDTH 
      ZCONSA = ZCONSA+ZDTHI5(JL)*(ZT0-PTSNM1M5(JL))*ZDTH
      ZCONSB = ZCONSB+ZDTHI5(JL)*(ZT0-PTSNM1M5(JL))*ZDTH
      PTSNM1M(JL) = PTSNM1M(JL)-ZDTHI5(JL)*(ZCONSA5(JL)+ZCONSB5(JL))*ZDTH

      ZHFLUX = ZHFLUX-ZDTHI5(JL)*ZDTHI5(JL)*ZDTHI
      ZHFLUXRF1 = ZHFLUXRF1-ZDTHI5(JL)*ZDTHI5(JL)*ZDTHI
      ZRS = ZRS+(ZT0-ZTSURF5(JL,1))*ZDTHI5(JL)*ZDTHI5(JL)*ZDTHI
      !*ZTSURF(JL) = ZTSURF(JL)-ZRS5(JL)*ZDTHI5(JL)*ZDTHI5(JL)*ZDTHI
      IF(YDSOIL%LESNICE)THEN
        PTIAM1M(JL,1) = PTIAM1M(JL,1)-ZRS5(JL)*ZDTHI5(JL)*ZDTHI5(JL)*ZDTHI
      ELSE
        PTSAM1M(JL,1) = PTSAM1M(JL,1)-ZRS5(JL)*ZDTHI5(JL)*ZDTHI5(JL)*ZDTHI
      ENDIF
    ENDIF
! SOLVE IMPLICIT TEMPERATURE EQUATION
    ZTSTARI = ZTSTARI+(PTSNM1M5(JL)+ZCONS25(JL) &
     & *(ZHFLUX5(JL)+ZHFLUXRF15(JL)+ZTSURF5(JL,1)*ZRS5(JL)))*ZTSTAR
    PTSNM1M(JL) = PTSNM1M(JL)+ZTSTARI5(JL)*ZTSTAR
    ZCONS2 = ZCONS2+ZTSTARI5(JL) &
     & *(ZHFLUX5(JL)+ZHFLUXRF15(JL)+ZTSURF5(JL,1)*ZRS5(JL))*ZTSTAR
    ZHFLUX = ZHFLUX+ZTSTARI5(JL)*ZCONS25(JL)*ZTSTAR
    ZHFLUXRF1 = ZHFLUXRF1+ZTSTARI5(JL)*ZCONS25(JL)*ZTSTAR
    ZRS = ZRS+ZTSTARI5(JL)*ZCONS25(JL)*ZTSURF5(JL,1)*ZTSTAR
    !*ZTSURF(JL) = ZTSURF(JL)+ZTSTARI5(JL)*ZCONS25(JL)*ZRS5(JL)*ZTSTAR
    IF(YDSOIL%LESNICE)THEN
      PTIAM1M(JL,1) = PTIAM1M(JL,1)+ZTSTARI5(JL)*ZCONS25(JL)*ZRS5(JL)*ZTSTAR
    ELSE
      PTSAM1M(JL,1) = PTSAM1M(JL,1)+ZTSTARI5(JL)*ZCONS25(JL)*ZRS5(JL)*ZTSTAR
    ENDIF

    ZRS = ZRS-ZCONS25(JL)*ZTSTARI5(JL)*ZTSTARI5(JL)*ZTSTARI
    ZCONS2 = ZCONS2-ZRS5(JL)*ZTSTARI5(JL)*ZTSTARI5(JL)*ZTSTARI
    !mass intercepted - recalculate from ZHFLUXRF1
    ZWFLUXRF = ZWFLUXRF+ZSSTAR*PTMST
    ZHFLUXRF1 = ZHFLUXRF1+ZWFLUXRF/RLMLT
    ZWFLUXRF = 0.0_JPRB

!  ACCOUNT FOR THROUGHFALL INTERCEPTION
    ! phase changes
    IF (ZHFLUXRF25(JL) < 0.0_JPRB) THEN
      ZHFLUXRF1  = 0.0_JPRB
    ENDIF
    IF (ZCOND55(JL) < ZCOND65(JL)) THEN
      ZCOND5 = ZCOND5+ZHFLUXRF1
    ELSE
      ZCOND6 = ZCOND6+ZHFLUXRF1
    ENDIF
    ZWFLUXRF = ZWFLUXRF+RLMLT*ZTMST*ZCOND6
    ZCONSA = ZCONSA+ZTMST*(ZT0-PTSNM1M5(JL))*ZCOND5
    PTSNM1M(JL) = PTSNM1M(JL)-ZTMST*ZCONSA5(JL)*ZCOND5
    ZRS = ZRS+ZT0*ZCOND5
    IF (ZCOND35(JL) < ZCOND45(JL)) THEN
      ZCOND3 = ZCOND3+ZWFLUXRF
    ELSE
      ZCOND4 = ZCOND4+ZWFLUXRF
    ENDIF
    ZWFLUXRF = 0.0_JPRB

    ZLWC(JL) = ZLWC(JL)+(1.0_JPRB-ZFUNC5(JL))*ZCOND4
    ZFUNC = ZFUNC-ZLWC5(JL)*ZCOND4
    ZWFLUXRF = ZWFLUXRF+PTMST*ZCOND3
    ! mass advected
    ZTSFC(JL) = ZTSFC(JL)+(PFRTI(JL,5)+PFRTI(JL,7))*ZWFLUXRF
    ZTSFL(JL) = ZTSFL(JL)+(PFRTI(JL,5)+PFRTI(JL,7))*ZWFLUXRF

    ! RESISTANCE FOR HEAT FLUX BETWEEN SNOW AND SOIL LAYER
    ZSNRES = ZSNRES-ZRS5(JL)*ZRS5(JL)*ZRS
    ZDSN = ZDSN+0.5_JPRB*ZSNRES/(RLAMICE*(PRSNM1M(JL)*ZHOICE)**RALAMSN)
! REAL SNOW DEPTH
    IF (ZCOND25(JL) < ZRDSNRESMAX) THEN
      ZCOND2 = ZCOND2+ZDSN
    ENDIF
    PSSNM1M(JL) = PSSNM1M(JL)+ZCOND2/(PRSNM1M(JL)*ZFRSN(JL))

! TERMS FOR IMPLICIT TEMPERATURE EQUATION
    ZCONSA = ZCONSA-PTMST*ZCONS2/(ZCONSA5(JL)+ZCONSB5(JL))**2
    ZCONSB = ZCONSB-PTMST*ZCONS2/(ZCONSA5(JL)+ZCONSB5(JL))**2
    IF (ZCOND15(JL) < ZCONSAMAX(JL)) THEN
      ZCOND1 = ZCOND1+ZCONSA
    ENDIF
    PSSNM1M (JL) = PSSNM1M (JL)+RHOCI*ZHOICE*ZIFRSN(JL)*ZCOND1
    ZDFUNC = ZDFUNC+RLMLT*ZLWC5(JL)*ZCONSB
    ZLWC(JL) = ZLWC(JL)+RLMLT*ZDFUNC5(JL)*ZCONSB

! SNOW LIQUID WATER - IMPLICIT
    ! ANALYTICAL FUNCTIONS - SNOW LIQUID WATER
    IF (PTSNM1M5(JL) < (ZT0-0.5_JPRB*RTEMPAMP) ) THEN
      ZDFUNC = 0.0_JPRB
      ZFUNC = 0.0_JPRB
    ELSEIF (PTSNM1M5(JL) < ZT0) THEN
      IF (LDREGSF) THEN
        ZDFUNC  = ZDFUNC/100.0_JPRB
      ENDIF
      PTSNM1M(JL) = PTSNM1M(JL)-RPI*RPI*ZITEMPAMP*ZITEMPAMP &
       & *SIN(RPI*(PTSNM1M5(JL)-ZT0)*ZITEMPAMP)*ZDFUNC
      PTSNM1M(JL) = PTSNM1M(JL)+RPI*ZITEMPAMP &
       & *COS(RPI*(PTSNM1M5(JL)-ZT0)*ZITEMPAMP)*ZFUNC
    ENDIF

    ! LWC IS RESCALLED CONSIDERING EFFECTIVE DEPTH (ZRDSNMAX)
    ! ZLWC(JL)=ZLWC(JL)*MIN(ZDSNR(JL),ZRDSNMAX)/ZDSNR(JL)
    ZDSNRP(JL) = ZDSNRP(JL)+ZLWC15(JL)*ZLWC(JL)/ZDSNR5(JL)
    ZDSNR(JL) = ZDSNR(JL)-ZLWC15(JL)*ZDSNRP5(JL)*ZLWC(JL)/(ZDSNR5(JL)**2)
    ZLWC(JL) = ZDSNRP5(JL)*ZLWC(JL)/ZDSNR5(JL)
    IF (ZDSNR5(JL) < ZRDSNMAX) THEN
      ZDSNR(JL) = ZDSNR(JL)+ZDSNRP(JL)
    ENDIF
    ZDSNRP(JL) = 0.0_JPRB

    ! SNOW LIQUID WATER CAPACITY
    PSSNM1M(JL) = PSSNM1M(JL) &
     & +(RLWCSWEB+(RLWCSWEC-RLWCSWEB)*ZCONSLWC(JL))*ZLWC(JL)
    ZLWC(JL) = 0.0_JPRB
    ! SNOW DEPTH (REAL)
    PSSNM1M(JL) = PSSNM1M(JL)+ ZDSNR(JL)/(PRSNM1M(JL)*ZFRSN(JL))
    ZDSNR(JL) = 0.0_JPRB

! PRELIMINARY SNOW MASS
   ! CORRECT FOR NEGATIVE VALUES OF SNOW MASS AFTER P-E
    IF (ZSSTARC5(JL) < 0.0_JPRB) THEN
      ZSSTAR  = 0.0_JPRB
    ENDIF
    PSSNM1M(JL) = PSSNM1M(JL)+ZSSTAR
    ZSSFC(JL) = ZSSFC(JL)+PTMST*ZSSTAR
    ZSSFL(JL) = ZSSFL(JL)+PTMST*ZSSTAR
    PEVAPTI(JL,5) = PEVAPTI(JL,5)+PTMST*PFRTI(JL,5)*ZSSTAR
    PEVAPSNW(JL) = PEVAPSNW(JL)+PTMST*PFRTI(JL,7)*ZSSTAR

!           NET HEAT FLUX AT SNOW SURFACE; THIS IS MULTIPLIED
!           BY FRACTION BECAUSE EQUATIONS APPLY TO TOTAL SNOW MASS
!           IN THE GRID SQUARE

    ZHFLUXP = ZHFLUXP+ZIFRSN(JL)*ZHFLUX

    PAHFSTI(JL,5) = PAHFSTI(JL,5)+PFRTI(JL,5)*ZHFLUXP
    PEVAPTI(JL,5) = PEVAPTI(JL,5)+PFRTI(JL,5)*RLSTT*ZHFLUXP
    PAHFSTI(JL,7) = PAHFSTI(JL,7)+PFRTI(JL,7)*ZHFLUXP
    PEVAPTI(JL,7) = PEVAPTI(JL,7)+PFRTI(JL,7)*RLVTT*ZHFLUXP
    PEVAPSNW(JL) = PEVAPSNW(JL)-PFRTI(JL,7)*RLVTT*ZHFLUXP
    PEVAPSNW(JL) = PEVAPSNW(JL)+PFRTI(JL,7)*RLSTT*ZHFLUXP
    PSSRFLTI(JL,5) = PSSRFLTI(JL,5)+PFRTI(JL,5)*ZHFLUXP
    PSSRFLTI(JL,7) = PSSRFLTI(JL,7)+PFRTI(JL,7)*ZHFLUXP
    PSLRFL(JL) = PSLRFL(JL)+(PFRTI(JL,5)+PFRTI(JL,7))*ZHFLUXP

    !*IF (LDLAND(JL)) THEN
    !*    PTSAM1M(JL,1) = PTSAM1M(JL,1) + ZTSURF(JL)
    !*    ZTSURF(JL)=0._JPRB
    !*ELSE
    !*  IF(YDSOIL%LESNICE)THEN
    !*    PTIAM1M(JL,1)=PTIAM1M(JL,1) + ZTSURF(JL)
    !*    ZTSURF(JL)=0._JPRB
    !*  ELSE
    !*    PTSAM1M(JL,1) = PTSAM1M(JL,1) + ZTSURF(JL)
    !*    ZTSURF(JL)=0._JPRB
    !*  ENDIF
    !*ENDIF

  ENDIF
ENDDO

!*         0.2. NEW SNOW T AND MASS INCLUDING GROUND HEAT FLUX AND MELTING.
!               -----------------------------------------------------------

DO JL=KIDIA,KFDIA
  PTSFL(JL) = PTSFL(JL)+ZGRIDFRAC(JL)*ZTSFL(JL)
  ZTSFL(JL) = 0.0_JPRB
  PTSFC(JL) = PTSFC(JL)+ZGRIDFRAC(JL)*ZTSFC(JL)
  ZTSFC(JL) = 0.0_JPRB
  PSSFL(JL) = PSSFL(JL)+ZGRIDFRAC(JL)*ZSSFL(JL)
  ZSSFL(JL) = 0.0_JPRB
  PSSFC(JL) = PSSFC(JL)+ZGRIDFRAC(JL)*ZSSFC(JL)
  ZSSFC(JL) = 0.0_JPRB
ENDDO

END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('SRFSN_LWINPSAD_MOD:SRFSN_LWIMPSAD',1,ZHOOK_HANDLE)

END SUBROUTINE SRFSN_LWIMPSAD
END MODULE SRFSN_LWIMPSAD_MOD
